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ABSTRACT 

A simple  and  relatively  efficient  method  for  simulating 
one-dimensional  and  two-dimensional  nonhomogeneous  Poisson  pro- 
cesses is  presented.  The  method  is  applicable  for  any  rate  function 
and  is  based  on  controlled  deletion  of  points  in  a Poisson  process 
whose  rate  function  dominates  the  given  rate  function.  In  its 
simplest  implementation,  the  method  obviates  the  need  for 
numerical  integration  of  the  rate  function,  for  ordering  of 
points,  and  for  generation  of  Poisson  variates. 


KEYWORDS:  Poisson  processes,  nonhomogeneous  Poisson  processes, 
one-dimensional  nonhomogeneous  Poisson  processes, 
two-dimensional  nonhomogeneous  Poisson  processes, 
simulation  of  point  processes,  thinning,  search 
and  detection,  acceptance-rejection,  order 
statistics. 


1 . INTRODUCTION 

The  one -dimensional  nonhomogeneous  Poisson  process  (see 
e.g.,  Cox  and  Lewis,  1966,  pp.  28-29;  £inlar,  1975,  pp.  94-101) 
has  the  characteristic  properties  that  the  numbers  of  points  in 
any  finite  set  of  nonoverlapping  intervals  are  mutually  independent 
random  variables,  and  that  the  number  of  points  in  any  interval 
has  a Poisson  distribution.  The  most  general  nonhomogeneous 
Poisson  process  can  be  defined  in  terms  of  a monotone  nondecreasing 
right -continuous  function  A(x)  which  is  bounded  in  any  finite 
interval.  Then  the  number  of  points  in  any  finite  interval,  for 
example  (0,x0J,  has  a Poisson  distribution  with  parameter 
* A(x0)  - A (0) . In  this  paper  it  is  assumed  that  A(x)  is 
continuous,  but  not  necessarily  absolutely  continuous.  The 
right  derivative  A(x)  of  A(x)  is  called  the  rate  function  of  the 
process;  A(x)  is  called  the  integrated  rate  function  and  has  the 
interpretation  that  for  x >_  0,  A(x)  - A ( 0 ) * E[N(x)],  where  N(x) 
is  the  total  number  of  points  in  (0,xj.  Note  that  \ (x)  may 
jump  at  points  at  which  A(x)  is  not  absolutely  continuous.  In 
contrast  to  the  homogeneous  Poisson  process,  i.e.,  A(x)  a constant 
(usually  denoted  by  A) , the  intervals  between  the  points  in  a 
one-dimensional  nonhomogeneous  Poisson  process  are  neither  inde- 
pendent nor  identically  distributed. 

Applications  of  the  one-dimensional  nonhomogeneous  Poisson 
process  include  modelling  of  the  incidence  of  coal-mining  disasters 
(Cox  and  Lewis,  1966),  the  arrivals  at  an  intensive 
care  unit  (Lewis,  1972),  transaction  processing  in  a data  base 
management  system  (Lewis  and  Shedler,  1976b) , occurrences  of  major 
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freezes  in  Lake  Constance  (Steinijans,  1976) , and  geomagnetic 
reversal  data  (Reyment,  1976) . The  statistical  analysis  of  trends 
in  a one -dimensional  nonhomogeneous  Poisson  process,  based  on  the 
assumption  of  an  exponential  polynomial  rate  function,  is  discussed 
by  Cox  and  Lewis  (1966)  ,-  Cox  (1972)  , Lewis  (1972)  , and  Lewis  and 
Shedler  (1976b). 

There  are  a number  of  methods  for  simulating  nonhomogeneous 
Poisson  process  which  we  review  briefly. 

(i)  Time-scale  transformation  of  a homogeneous  (rate  one) 

Poisson  process  via  the  inverse  of  the  (continuous)  integrated 
rate  function  A(x)  constitutes  a general  method  for  generation 
of  the  nonhomogeneous  Poisson  process  (cf.,  £inlar,  1975 
pp.  96-97).  This  method  is  based  on  the  result  that  X. ,X2,..,, 
are  the  points  in  a nonhomogeneous  Poisson  process  with  continuous 
integrated  rate  function  A (x)  if  and  only  if  XJ  * A (x^ , 

X2  * A(X2),  ...  , are  the  points  in  a homogeneous  Poisson  process 
of  rate  one.  The  time-scale  transformation  method  is  a direct 
analogue  of  the  inverse  probability  integral  transformation 
method  for  generating  (continuous)  nonuniform  random  numbers. 

For  many  rate  functions,  inversion  of  A(x)  is  not  simple 
and  must  be  done  numerically;  cf.,  Gilchrist  (1977)  and  Patrow 
(1977) . The  resulting  algorithm  for  generation  of  the  non- 
homogeneous Poisson  process  may  be  far  less  efficient  than 
generation  based  on  other  methods;  see  e.g.,  Lewis  and  Shedler 
(1976a,  1977)  and  Patrow  (1977)  for  discussion  of  special 
methods  for  efficiently  generating  the  nonhomogeneous  Poisson 
process  with  log-linear  and  log-quadratic  rate  functions. 
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(ii)  A second  general  method  for  generating  a nonhomogeneous  Poisson 
process  with  integrated  rate  function  A (x)  is  to  generate  the 
intervals  between  points  individually,  an  approach  which  may 
seem  more  natural  in  the  event-scheduling  approach  to  simulation. 
Thus,  given  the  points  X^  = x^,  X2  * x2,...,  X^  = x^,  with 
Xi  < X2  <•••<  X.,  the  interval  to  the  next  point,  X^+1  - X^ 
is  independent  of  x^,  ...  , x^^  and  has  distribution  function 

F(x)  = 1 - expt-tAfx^  + x)  - A(x^))J.  It  is  possible  to  find  the 
inverse  distribution  function  usually  numerically, 

and  generate  X^+1  - X^  according  to  - X^  = F**1  (U^)  , where 

is  a uniform  random  number  on  the  interval  (0,1) . Note, 
however,  that  this  not  only  involves  computing  the  inverse  distri- 
bution function  for  each  interval  X^+^-  X^,  but  that  each  distribution 
has  different  parameters  and  possibly  a different  form.  An 
additional  complication  is  that  X^+1  - X^  is  not  necessarily  a 
proper  random  variable,  i.e.,  there  may  be  positive  probability 
. that  Xi+1  - xi  is  infinite.  It  is  necessary  to  take  this  into 
account  for  each  interval  x^+1”  X^  before  the  inverse  probability 
integral  transformation  is  applied.  The  method  is  therefore  very 
inefficient  with  respect  to  speed,  more  so  than  the  time-scale 
transformation  method. 

(iii)  In  a third  method, simulation  of  a nonhomogeneous  Poisson 

process  in  a fixed  interval  (0,xQ]  can  be  reduced  to  the 
generation  of  a Poisson  number  of  order  statistics  from  a fixed 
density  function  by  the  following  result  (cf.,  Cox  and  Lewis, 

1966,  p.  45).  If  X^ , X2,  ...  , XR  are  the  points  of  the 
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nonhomogeneous  Poisson  process  in  (0,xQ],  and  if  N(xQ)  = n, 

then  conditional  on  having  observed  n(>  0)  points  in  (0,Xg], 

the  are  distributed  as  the  order  statistics  from  a sample 

of  size  n from  the  distribution  function  (A(x)  - A(0)  }/{A(xQ)  -A(0)}, 

defined  for  0 < x £ Xq.  Generation  of  the  nonhomogeneous  Poisson 

process  based  on  order  statistics  is  in  general  more  efficient 

(with  respect  to  speed)  than  either  of  the  previous  two  methods. 

Of  course,  a price  is  paid  for  this  greater  efficiency.  First, 

it  is  necessary  to  be  able  to  generate  Poisson  variates,  and 

second,  more  memory  is  needed  than  in  the  interval -by-interval 

method  in  order  to  store  the  sequence  of  points.  Enough  memory 

must  be  provided  so  that  with  very  high  probability  the  random 

number  of  points  generated  in  the  interval  can  be  stored.  Recall 

that  the  number  of  points  in  the  interval  (0,xQ]  has  a Poisson 

distribution  with  mean  « A(xQ)  ••  A (0)  . Memory  of  size,  e.g.  , 

1/2 

M0  + 4Wq  will  ensure  that  overflow  will  occur  on  the  average 
in  only  one  out  of  approximately  every  40,000  realizations.  This 
probability  is  small  enough  so  that  in  the  case  of  overflow,  the 
realization  of  the  process  can  generally  be  discarded. 

(iv)  Again,  there  is  a very  particular  and  very  efficient 

method  for  simulation  of  nonhomogeneous  Poisson  processes  with 
log-linear  rate  function  (Lewis  and  Shedler,  1976a)  which,  at 
the  cost  of  programming  complexity  and  memory,  can  be  used  to 
obtain  an  efficient  simulation  method  for  other  rate  functior  ,, 
as  in  Lewis  and  Shedler  (1977) . 
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In  this  paper  a new  method  is  given  for  simulating  a non- 
homogeneous  Poisson  process  which  is  not  only  conceptually  simple, 
but  is  also  computationally  simple  and  relatively  efficient.  In 
fact,  at  the  cost  of  some  efficiency,  the  method  can  be  applied 
to  simulate  the  given  nonhomogeneous  Poisson  process  without  the 
need  for  numerical  integration  or  routines  for  generating  Poisson 
variates . Used  in  conjunction  with  the  special  methods  given  by 
Lewis  and  Shedler  (1976a,  1977),  the  method  can  be  used  to  generate 
quite  efficiently  nonhomogeneous  Poisson  processes  with  rather 
complex  rate  function,  in  particular  combinations  of  long-term 
trends  and  fixed-cycle  effects.  The  method  is  also  easily  extended 
to  the  problem  of  generating  the  two-dimensional  nonhomogeneous 
Poisson  process. 


2.  SIMULATION  OF  ONE -DIMENSIONAL  NONHOMOGENEOUS  POISSON  PROCESSES 

Simulation  of  a nonhomogeneous  Poisson  process  with  general 

rate  function  A(x)  in  a fixed  interval  can  be  based  on  thinning 

* 

of  a nonhomogeneous  Poisson  process  with  rate  function  A (x)  A(x). 

The  basic  result  is 


THEOREM  1.  Consider  a one -dimensional  nonhomogeneous  Poisson 

* * 

process  {N  (x)  ?x  _>  0}  with  rate  function  A (x) , so  that  the 

* 

number  of  points,  N (xQ) , in  a fixed  interval  (0,xQ]  has  a 
Poisson  distribution  with  parameter  Uq  » A*(xQ)  - A*(0).  Let 


* 

<1 


X,  , X. 
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Xjj  *^x  j be  the  points  of  the  process  in  the  interval 


5 


( 0 , Xg ] . Suppose  that  for  0 < x < xQ,  X(x)  < X*(x).  For 

i = l,2,...,n,  delete  the  point  X?  with  probability 

* * * 

1 - X(X.)/X  (Xi) ; then  the  remaining  points  form  a nonhomogeneous 
Poisson  process  {N(x):x  0}  with  rate  function  X(x)  in  the 

interval  (0,xQ]. 


Proof , Since  {N  (x)  jx  0)  is  a nonhomogeneous  Poisson  process, 
and  points  are  deleted  independently,  it  is  clear  that  the  number 
of  points  in  {N(x);x  0}  in  any  set  of  nonoverlapping  intervals 

are  mutually  independent  random  variables.  Thus,  it  is  sufficient 
to  show  that  the  number  of  points  N(a,b)  in  (N(x)  ;x  0} 
in  an  arbitrary  interval  (a,b]  with  0 < a < b £ xQ  has  a 
Poisson  distribution  with  parameter  A(b)  - A (a). 

Observe  that  with  p(a,b)  « (A(b)  - A (a) }/{ A* (b)  - A*(a)}, 
we  have  the  conditional  probability 


P{N (a ,b) 


n|N*(a,b)*k) 


1 

{p(a,b)}n{l-p(a,b)}k”n 


0 


if 

if 

and 

if 

and 


n * k = 0 

k n > 0 
k > 1 

n > 1 
k < n . 


(1) 


Equation  (1)  is  a consequence  of  the  well-known  result  that, 
conditional  on  n ( >0 ) points  in  the  interval  (a,bj,  the  joint 

it 

density  of  the  n points  in  the  process  {N  (x) :x  > 0}  is 
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• • # 


A*(x1) 


A*(x  )/{A*(b)  - A*(a))n.  The  desired  result  is 
n 


obtained  in  a straightforward  manner  from  Equation  (1)  by  removing 
the  condition. 


Theorem  1 is  the  basis  for  the  method  of  simulating  non- 
homogeneous  Poisson  processes  given  in  this  paper. 

ALGORITHM  l.  One-dimensional  nonhomogeneous  Poisson  process. 

if 

1.  Generate  points  in  the  nonhomogeneous  Poisson  process  (N  (x)}with 

h 

rate  function  A (x)  in  the  fixed  interval  (0,x»].  If 

* * 

the  number  of  points  generated,  n , is  such  that  n =0, 
exit;  there  are  no  points  in  the  process  {N(x)}. 

Hr  * * 

2.  Denote  the  (ordered)  points  by  X^,  X2,  ...  , 3^*.  Set 
i « 1 and  k ■ 0. 

3.  Generate  U.  , uniformly  distributed  between  0 and  1.  If 

* * * * 

U.  £ A(X^)/A  (X.),  set  k equal  to  k+1  and  X^  » X^. 

* 

4.  Set  i equal  to  i+1.  If  i £ n , go  to  3. 

5.  Return  X^ , X2,  ...  , Xn,  where  n = k,  and  also  n. 

In  the  case  where 

* * * 

(i)  {N  (x) } is  a homogeneous  Poisson  process  with  A (x)  * A ; 

(ii)  the  minimum  of  A(x),  say  £,  is  known,  and 

(iii)  generation  of  uniformly  distributed  variates  is  computa- 
tionally costly, 

* 

considerable  speedup  can  be  obtained  by  noting  that  Xi  is  always 

* 

accepted  if  IL  £ A/ A . This  obviates,  in  some  cases,  computation 
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of  x(x) , which  is  the  main  source  of  inefficiency  in  the  algorithm. 

h 

Moreover,  in  this  case  X U^/X  can  be  used  as  the  next  uniformly 
distributed  variate. 

The  method  of  thinning  in  this  simple  form,  i.e., 

X* (x)  « X*  > max.  , „ , „ X(x),  can  also  be  used  to  provide  an 

*"  U N X N Xq 

algorithm  for  generating  a nonhomogeneous  Poisson  process  on  an 
interval -by-interval  basis,  as  discussed  in  subsection  (ii)  of 
Section  1.  The  interval  to  the  next  point  - X^  is 

obtained  by  generating  and  cumulating  exponential  (a*)  random 


numbers  E 


1 * ' * * * ' 


until  for  the  first  time 
* * * 

Uj  £ X^  + E^  + •••  + Ej)/X  , where  the  Uj  are  independent 
uniform  random  numbers  between  0 and  1.  This  algorithm  is 
considerably  simpler  than  the  interval-by-interval  algorithm 
of  Section  1 since  it  requires  no  numerical  integration, 
only  the  availability  of  uniform  random  numbers. 
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3.  DISCUSSION  OF  THE  METHOD  OF  THINNING 

(i)  Relationship  to  acceptance-rejection  method 

The  method  of  thinning  of  Algorithm  1 is  essentially  the 
obverse  of  the  conditional  method  of  Section  1,  using  condition- 
ing and  acceptance-rejection  techniques  to  generate  the  random 
variables  with  density  function  X(x)/{A(x)  - A ( 0 ) } (Lewis 
and  Shedler,  1977,  Algoritnm  3).  The  differences  are  subtle, 
but  computationally  important.  In  the  acceptance-rejection 
metnod,  it  is  first  necessary  to  generate  a Poisson  variate 
with  mean  * A(xQ)  - A(0),  and  this  involves  an  integration 

of  the  rate  function  X (x) . Then  the  Poisson  (y0)  number,  n, 
of  variates  generated  by  acceptance-rejection  must  be  ordered  to 
give  X^,  X2,  ...  , Xn. 

(ii)  Simplest  form  of  the  thinning  algorithm 

In  the  simplest  form  of  the  method  of  thinning,  X*(x)  is 
taken  to  be  a constant  X , so  that,  for  instance,  the  Doints 
Xl'  X2'  ' Xn*  can  be  generated  by  cumulating  exponential  (X  ) 

variates  until  the  sum  is  greater  than  xQ  (cf.,  Lewis  and 
Shedler,  1976a,  Algorithm  1) . Thinning  is  then  applied  to  the 
generated  points.  No  ordering,  no  integration  of  X(x)  and  no 

generator  of  Poisson  variates  is  required.  Of  course  for  both 
algorithms  to  be  efficient,  computation  of  X (x)  and  X*(x) 
must  be  easy  relative  to  computation  of  the  inverse  of  A (x) . 
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(iii) 


Efficiency 

For  the  thinning  algorithm  (as  well  as  the  algorithm  based 

on  conditioning  and  acceptance-rejection)  efficiency,  as  measured 

* 

by  the  number  of  points  deleted,  is  proportional  to  viq/Mq  = 

{A(x0)  - A (0 ) }/{A*  (xQ)  - A*  (0) };  this  is  the  ratio  of  the 

* * 

areas  between  0 and  Xq  under  A(x)  and  A (x) . Thus,  A (x) 
should  be  as  close  as  possible  to  A(x)  consistent  with  ease 
of  generating  the  nonhomogeneous  Poisson  process  {N*(x):x  >_  0}. 

(iv)  An  example:  fixed  cycle  plus  trend 

To  illustrate  the  applicability  of  the  thinning 
algorithm,  consider  its  use  in  conjunction  with  the  algorithms 
given  by  Lewis  and  Shedler  (1976a,  1977)  for  log-linear  and 
log-quadratic  rate  functions.  Assume  that  it  is  necessary  to 
simulate  a nonhomogeneous  Poisson  process  whose  rate  function 
increases  quadratically  with  time  but  also  has  a fixed-period 
cycle,  e.g., 


2 

A(x)  » exp{aQ  + a^x  + oijX  + K sin  (u)Qx  * 9)}, 

0 < x < xQ;  K > 0,  0 < 0 £ 2tt  , > 0 . 

This  is  the  model  found  by  Lewis  (1972)  for  arrivals  at  an 
intensive  care  unit,  where  there  is  a strong  time-of-day  effect. 

Thus  if  u)Q  = 2tt/Tq,  then  the  period  Tq  * 1 day.  Computation 

—l  * 

of  A (♦)  is  difficult.  To  determine  A (x) , note  that 

* 2 

A(x)  £ A (x)  = exp{aQ  + K + a^x  + a^x  } , 

10 


and  therefore 


A(x)/A*(x)  = exp[K{l  - sin(tt)Qx  +0)}]  . 

Thus  in  step  3 of  Algorithm  1,  tL  is  compared  to 

it 

exp[K(l  - sindgX^  + 6)}].  Equivalently,  if  unit  exponential 

variates  E^  are  available,  it  is  faster  to  compare  E^  to 

K{1  - sin(u0X*  + 0)},  accepting  X^  if  E^  > K{1  - sin(w0X*  + 0)}. 

The  main  computational  expense  here  is  generation  of  the  Ei 

* 

and  computation  of  the  sine  function,  both  n times.  The  expense 

involved  in  computation  of  the  sine  function  can  be  reduced  by 

* 

noting  that  the  point  X^  is  always  accepted  if  E^  is  greater 
than  2K.  This  will  be  a great  saving  if  the  cyclic  effect  is 
minor  (K  small).  The  number  of  E^  generated  can  be  reduced 

by  noting  that  if,  in  one  step  of  the  algorithm,  Ei  is  observed 

* 

to  be  greater  than  6,  then  E^  « - 6 can  be  used  as  an 

(independent)  unit  exponential  variate  in  the  next  step.  The 
above  procedure  can  be  extended  to  the  case  of  a trend  with 
two  fixed-period  cycles,  e.g.,  a time-of-day  and  a time-of-week 
effect. 


4.  SIMULATION  OF  TWO-DIMENSIONAL  HOMOGENEOUS  POISSON  PROCESSES 

The  two-dimensional  homogeneous  Poisson  process  (of  rate  X > 0) 
is  defined  by  the  properties  that  the  numbers  of  points  in  any  finite 
set  of  nonoverlapping  regions  having  areas  in  the  usual  geometric 
sense  are  mutually  independent,  and  that  the  number  of  points  in  any 
region  of  area  A has  a Poisson  distribution  with  mean  AA; 
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see,  e.g.,  Karlin  and  Taylor  (1975),  pp.  31-32.  Note  that  the 
number  of  points  in  a region  R depends  on  its  area,  but  not  on  its 
shape  or  location.  The  homogeneous  Poisson  process  arises  as  a 
limiting  two-dimensional  point  process  with  respect  to  a number 
of  limiting  operations;  cf.,  Goldman  (1967a, b) . Properties  of 
the  process  are  given  by  Miles  (1970)  . Applications  of  the  two- 
dimensional  homogeneous  Poisson  process  to  problems  in  ecology 
and  forestry  have  been  discussed  by  Thompson  (1955)  and  Holgate 
(1972)  . The  model  also  arises  in  connection  with  naval  search 
and  detection  problems. 

In  considering  the  two-dimensional  homogeneous  Poisson 
process,  projection  properties  of  the  process  depend  quite 
critically  on  the  geometry  of  the  regions  considered.  These 
projection  properties  are  simple  for  rectangular  and  circular 
regions,  and  make  simulation  of  the  homogeneous  process  quite 
easy.  We  consider  these  two  cases  separately. 

(i)  Homogeneous  Poisson  Processes  in  a Rectangle 

The  following  two  theorems  form  the  basis  for  simulation 
of  the  two-dimensional  homogeneous  Poisson  process  in  a 
rectangle . 

THEOREM  2.  Consider  a two-dimensional  homogeneous  Poisson  process 
of  rate  A,  so  that  the  number  of  points  in  a fixed  rectangle 
r « { (x,y) : 0 < x < xQ,  0 < y £ yQ}  has  a Poisson  distribution 
with  parameter  AxQy0.  If  (X^Y^),  (X2,Y2)#  • ••  • denotcJ 

the  position  of  the  points  of  the  process  in  R,  labelled  so  that 
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form  a one-dimensional 


Xx  < X2  < •••  < XN,  then  Xr  X2<  ...  , Xjj 

homogeneous  Poisson  process  on  0 < x £ xQ  of  rate  XyQ.  If  the 

points  are  relabelled  (Xj,Y£),  (X£,Y2) , ...  / (X^Y^)  so  that 

Y1  < Y2  < < YN'  then  Yl'  Y2'  •••  ' YN  form  a one-dimensional 

homogeneous  Poisson  process  on  0 < y <_  of  rate  XxQ. 

Proof . The  number  of  points  in  an  interval  on  the  x-axis , say, 

(a,b]  is  the  number  of  points  in  the  rectangle  bounded  by  the 
lines  x ■ a,  x*b,  y=0,  and  y ■ yQ.  This  number  is  therefore 
independent  of  the  number  of  points  in  any  similar  nonoverlapping 
rectangle  bounded  on  the  x-axis  by  x * a',  x **  b',  i.e.,  the 
number  of  points  in  the  interval  (a',b'].  This  establishes  the 
independent  increment  property  for  a one-dimensional  Poisson 
process.  The  Poisson  distribution  of  the  number  of  points  in 
(a,b]  follows  from  the  fact  that  it  is  equal  to  the  number  of 
points  in  the  rectangle  bounded  by  x«a,  x*b,  y«0,  and 
y » y0,  and  the  latter  has  a Poisson  distribution  with  parameter 
AyQ(b-a).  An  analogous  argument  shows  that  the  process  formed 
on  the  y-axis  by  Y|,  Y£,  ...  , Y^  is  Poisson. 


Conditional  properties  of  the  Poisson  process  in  a rectangle 
are  established  next.  The  important  thing  to  note  is  that  while 
the  processes  obtained  by  projection  of  the  points  onto  the  x 
and  y axes  are  not  independent,  there  is  a type  of  conditional 
independence  which  makes  it  easy  to  simulate  the  two-dimensional 
process . 
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THEOREM  3.  Assume  that  a two-dimensional  homogeneous  Poisson 
process  of  rate  X is  observed  in  a fixed  rectangle 
R = Hx,y)  • 0 < x£Xq,  0 < y <_  Yq)  * so  that  the  number  of 
points  in  R,  N(R),  has  a Poisson  distribution  with  parameter 
Xx0y0.  If  N(R)  ® n > 0 and  if  (X1,Y1) , <X2 ,Y2) , . . . , (Xn,Yn) 

denote  the  points,  labelled  so  that  X^  < X2  < •••  < Xn,  then  con- 
ditional on  having  observed  n points  in  R,  the  X^  X2,...,Xn 

are  uniform  order  statistics  on  0 < x < x„,  and  Y_ , Y_ , ...  , Y 

— U 1 2 n 

are  independent  and  uniformly  distributed  on  0 < y £ y , 
independent  of  the  X^. 

Proof . If  there  are  N points  in  the  rectangle,  form  N 
vertical  strips,  from  0 to  yQ  and  from  Xi  to 
such  that  each  strip  contains  only  one  of  the  N points.  The 
position  of  Y^  on  the  vertical  line  through  Xi  is  that  of 
an  event  in  a Poisson  process  of  rate  Xdx^,  given  that  only 
one  event  occurs.  But  this  means  that  Y^  is  uniformly  dis- 
tributed between  0 and  y^ . Moreover,  this  is  true  irrespective 
of  where  X^  occurs;  therefore  Y^  is  independent  of  xi . 

Also,  occurrences  in  all  N strips  are  independent,  and  therefore 
Y^  is  independent  of  the  other  Yj  and  X..  positions,  j + i. 

Tnus  the  Y^  are  a random  sample  of  size  N from  a uniform 
(0,yQ)  distribution,  independent  of  the  X^.  Now  condition 

on  N * n (>  0);  since  by  Theorem  2 the  X^  form  a Poisson  process 
they  are,  by  well-known  results,  order  statistics  from  a uniform 
(0,Xq)  sample  and  are  independent  of  the  fixed  size  Y^  popula- 
tion* .thiis  the  pairs  (X^,Y^)  are  mutually  independent. 
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COROLLARY.-  Denote  the  Poisson  points  by  (X^Y^,  (X2,Y2),  . . ., 


where  the  index  does  not  necessarily  denote  an  ordering  on  either 
axis.  Conditionally,  the  pairs  (X]/Yi) , . . . , (X^Yjj)  are  inde- 
pendent random  variables.  Furthermore,  for  each  pair  , X^  is 
distributed  uniformly  between  0 and  xQ,  independently  of  Y^ 

which  is  uniformly  distributed  between  0 and  yQ. 

• * 

From  the  two  theorems,  the  following  simulation  procedure 
is  obtained. 

ALGORITHM  2.  Two-dimensional  homogeneous  Poisson  process  in  a 
rectangle. 

1.  Generate  points  in  the  one-dimensional  homogeneous  Poisson 
process  of  rate  XyQ  on  (0,xQ].  If  the  number  of 
points  generated,  n,  is  such  that  n « 0,  exit;  there 

are  no  points  in  the  rectangle. 

2.  Denote  the  points  generated  by  X1  < X2  < •••  < XR. 

3.  Generate  Yj.,  Y2,  ...  , Y as  *ndePendent'  uniformly  dis- 
tributed random  numbers  on  (0,yQJ. 

4.  Return  (X^Y^,  (X2,Y2),  ...  , (Xn,Yn)  as  the  coordinates 
of  the  two-dimensional  homogeneous  Poisson  process  in  the 
rectangle,  and  n. 

Note  that  generation  of  the  points  , X2 , . . . , X^  in  Steps  1 
and  2 can  be  accomplished  by  cumulating  exponential (XyQ)  random 
numbers.  Alternatively,  after  generating  a Poisson  random  number 
N = n (with  parameter  XxQy0) , n independent,  uniformly  distributed 
random  numbers  on  (0,xQ]  can  be  ordered;  see  Lewis  and  Shedler 
(1976a),  p.  502. 
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Another  algorithm  for  generation  of  the  two-dimensional 
Poisson  process  in  a rectangle  can  be  based  on  the  Corollary 
to  Theorem  3. 

(ii)  Homogeneous  Poisson  Processes  in  a Circle 

The  following  two  theorems  form  the  basis  for  simulation 
of  tne  two-dimensional  homogeneous  Poisson  process  in  a fixed 
circle  of  radius  Tq. 

Fix  the  origin  and  initial  line  of  polar  coordinates  r 
and  6 so  that  the  origin  is  the  center  of  the  circle  and  the 
initial  line  is  horizontal.  We  consider  the  projection  of  the 
points  (R^,0^),  of  the  Poisson  process  circularly  onto  the 
r-axis  (R^)  and  radially  onto  the  circumferential  e-axis  (6^. 

t 

The  number  of  points  projected  onto  the  r-axis  in  the  interval 

(0,r],  where  r £ rQ,  is  the  number  of  points  in  the  circle  of 
~ 2 

radius  r and  area  irr  ; thus  the  number  of  points  in  (0,rj 

2 

has  a Poisson  distribution  with  parameter  Airr  . Consequently, 

if  the  projection  process  on  the  r-axis  is  a Poisson  process, it 

2 

must  have  integrated  rate  function  A(r)  = Airr  , with  A(0)  = 0. 

Similarly, the  number  of  points  on  the  circumferential  arc 
of  the  fixed  circle  (radius  Tq)  from  0 to  0 is  the  number 
of  points  in  the  sector  of  the  circle  defined  by  radial  lines 
at  angles  0 and  0;  thus, the  number  of  points  on  the  arc  from 

2 ft 

0 to  0 has  a Poisson  distribution  with  parameter  A7rrQ  x 
2 

* 0ArQ/2.  Accordingly,  if  the  projection  process  on  the  0-axis 
is  a Poisson  process,  it  must  have  integrated  rate  function 
A (6)  * OArg/2,  with  A (0)  = 0. 
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We  now  assert  that  the  projection  processes  are  in  fact 
Poisson  processes.  Since  proofs  of  these  theorems  are  directly 
analogous  to  the  proofs  of  Theorems  2 and  3,  they  are  omitted. 

THEOREM  4.  Consider  a two-dimensional  homogeneous  Poisson  process 

of  rate  X so  that  the  number  N of  points  in  a fixed  circular  area 

2 

C of  radius  r^  and  area  ttVq  has  a Poisson  distribution 

2 

with  parameter  X7rr0.  If  (R^^),  (R2,62),  ...  , (R^e^) 
denote  the  points  of  the  process  in  C,  labelled  so  that 

R1  < R2  < ***  < V tnen  R2'  ***  ' ^ a one-dimensional 

nonhomogeneous  Poisson  process  on  0 £ r £ Tq  with  rate  function 

X(r)  « 2irXr . If  the  points  are  relabelled  (R^,e2), 

. ..,  so  that  ej  < ej  < •••  < 0j|,f  then  6^, 

form  a one-dimensional  homogeneous  Poisson  process  on  0 < 9 £ 2tt 
2 

of  rate  XrQ/2. 

THEOREM  5.  Assume  that  a two-dimensional  Poisson  process  of 

rate  X is  observed  in  a fixed  circular  area  C of  radius  rQ 

so  that  the  number  of  points  in  C,  N(C),  has  a Poisson  distri- 

2 

bution  with  parameter  XirrQ.  If  N(C)  * n > 0 and  if 

^,0^,  (R2»92)»  * (Rn'en)  with  Ri  < V ”*  < Rn 

denote  the  points,  then  conditional  on  having  observed  n points  in 

C,  the  R^,  R2,  ...  , Rr  are  order  statistics  from  the  density 
2 

f(r)  * 2r/rQ  concentrated  on  0 £ r < rQ,  and  6^,  02,  ...  , 0R 
are  independent  and  uniformly  distributed  on  0 < e < 2tt, 
independent  of  the  R^ . 
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These  theorems  lead  to  the  following  simulation  procedure. 


ALGORITHM  3.  Two-dimensional  homogeneous  Poisson  process  in  a 
circular  area. 

1.  Generate  n as  a Poisson  random  number  with  parameter 

2 

Affrg.  If  n « 0,  exit;  there  are  no  points  in  C. 

2.  Generate  n independent  random  numbers  having  density 

function  f (r)  = 2r /ti.  and  order  to  obtain  R.  < R»  < • • • < R . 

u i 2 n 

3.  Generate  62,...  ' ®n  independent,  uniformly  distributed 

random  numbers  on  (0,2tt]. 

4.  Return  (R^e^),  (R2,62),  ...  , (Rn»6n)>  and  n. 

Note  that  the  wedge-shaped  density  2r/rg  can  be  generated  by 
scaling  the  maximum  of  two  independent  uniform  (0,1)  random 
numbers . 

Direct  generation  of  homogeneous  Poisson  points  in 
non-circular  or  non-rectangular  regions  is  difficult.  The 
processes  obtained  by  projection  of  the  points  on  the  two  axes 
are  nonhomogeneous  Poisson  processes  with  complex  rate  functions 
determined  by  the  geometry  of  the  region.  However,  the  conditional 
independence  which  is  found  in  circular  and  rectangular  regions 
(Theorems  3 and  5)  for  the  processes  on  the  two  axes  is  not 
present.  In  particular,  given  that  there  are  n points 
(X1,Y1) , . . . , (Xn,Yn)  in  a non-rectangular  region,  the  pairs  (X^Y^ 
are  mutually  independent,  but  XA  is  in  general  not  independent 
of  Yi,  i = 1, . . . ,n.  Therefore,  it  is  simpler  to  enclose  the 
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region  in  either  a circle  or  a rectangle,  generate  a homogeneous 
Poisson  process  in  the  enlarged  area,  and  subsequently  exclude 
points  outside  of  the  given  region. 


5.  SIMULATION  OF  TWO-DIMENSIONAL  NONHOMOGENEOUS  POISSON  PROCESSES 
The  two-dimensional  nonhomogeneous  Poisson  process 
{N(x,y):x  > 0,  y >_  0}  is  specified  by  a positive  rate  function 
X (x,y)  which  for  simplicity  is  assumed  here  to  be  continuous. 

Then  the  process  has  the  characteristic  properties  that  the  numbers 

of  points  in  any  finite  set  of  nonoverlapping  regions  having 
areas  in  the  usual  geometric  sense  are  mutually  independent,  and 
that  the  number  of  points  in  any  such  region  R has  a Poisson 
distribution  with  mean  A(R);  here  A (R)  denotes  the  integral 
of  A(x,y)  over  R,  i.e.,  over  the  entire  area  of  R. 

Applications  of  the  two-dimensional .nonhomogeneous  Poisson 
process  include  problems  in  forestry  and  naval  search  and 
detection.  The  use  of  the  process  as  a model  for  the  pattern  of 
access  to  the  storage  subsystem  of  a computer  system  will  be 
reported  elsewhere.  Detection  and  statistical  analysis  of  trends 
in  the  two-dimensional  nonhomogeneous  Poisson  process  is  discussed 
by  Rantschler  (1973) . 

Theorem  1 dealing  with  thinning  of  one-dimensional  non- 
homogeneous Poisson  processes  generalizes  to  two-dimensional 

nonhomogeneous  Poisson  processes.  Thus,  suppose  that 

* 

A(x,y)  £ A (x,y)  in  a fixed  rectangular  region  of  the  plane. 

If  a nonhomogeneous  Poisson  process  with  rate  function  A (x,y) 
is  thinned  according  to  A(x,y)/A*(x,y)  (i.e.,  each  point 
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I 

\ 

f 

Xi'Yi^  deleted  independently  if  a uniform  (0,1)  random 
number  Ui  is  greater  than  X (X^Y^/A*  (Xi,Yi) ) , the  result 
is  a nonhomogeneous  Poisson  process  with  rate  function  A(x,y). 

1 The  proof  is  a direct  analogue  of  the  proof  for  the  one-dimensional 

case. 

The  nonhomogeneous  Poisson  process  with  rate  function 
A(x,y)  in  an  arbitrary  but  fixed  region  R can  be  generated  by 

! 

enclosing  the  region  R either  in  a rectangle  or  a circle,  and 

applying  Algorithm  2 or  Algorithm  3.  The  following  procedure 

assumes  that  the  region  R has  been  enclosed  in  a rectangle  R* , 

* 

and  that  A a max{A(x,y):  x,y  * R)  has  been  determined;  here 

the  bounding  process  is  homogeneous  with  rate  A*  in  the  rectangle  R*. 

ALGORITHM  4.  Two-dimensional  nonhoroogenous  Poisson  process. 

3.  Using  Algorithm  2,  generate  points  in  the  homogeneous 

* * 

Poisson  process  of  rate  A in  the  rectangle  R . If 

* * 

the  number  of  points,  n , is  such  that  n « 0,  exit;  there 

are  no  points  in  the  nonhomogeneous  Poisson  process. 

* 

2.  From  the  n points  generated  in  1,  delete  the  points  that 
are  not  in  R,  and  denote  the  remaining  points  by 

<xI'*l)'  <X2'V <Xm'V  with  h < x2  < •"  < V 

Set  i ■ 1 and  k ■ 0. 

3.  Generate  uniformly  distributed  between  0 and  1.  If 

Oi  < A(X*,Y*)/A*,  set  k = k+1,  XR  = X*  and  Yk  = Y* . 

* 

4.  Set  i equal  to  i+1.  If  i £ m , go  to  3. 

5.  Return  (X-^Y^) , (X2,Y2),  ...  , (Xn/Yn)  , where  n = k,  and  n. 
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It  is  not  necessary  that  the  bounding  process  have  a constant  rate 
* 

A . Theroems  2 and  4 can  be  generalized  to  certain  cases  where 
the  process  is  nonhomogeneous  (cf.,  Bartlett,  1974) , for  instance 
A(x,y)  - p (x)  <p(y).  Thus,  a tighter  bounding  process  which  is 
nonhomogeneous  may  possibly  be  obtained.  It  is  not  simple  to 
see  how  much  efficiency  could  be  gained  by  doing  this,  as  opposed 
to  using  a two-dimensional  homogeneous  Poisson  process  for  the 
bounding  process.  Again,  as  in  the  one-dimensional  case,  savings 
in  computing  A(x,y)  can  be  obtained  by  computing  its  minimum 
beforehand,  and  the  U^'s  can  be  reused  by  scaling. 

6.  COMPARISONS  AND  CONCLUDING  REMARKS 

The  method  of  thinning  presented  in  this  paper  for  simulating 
one-dimensional  and  two-dimensional  nonhomogeneous  Poisson  processes 
with  given  rate  function  can  be  carried  out  in  a computationally 
simple  way  by  using  a bounding  process  which  is  homogeneous  with  a 
rate  function  equal  to  the  maximum  value  of  the  given  rate  function. 
No  numerical  integration,  ordering  or  generation  of  Poisson 
variates  is  required,  only  the  ability  to  evaluate  the  given  rate 
function.  The  thinning  algorithm  appears  to  be  particularly 
attractive  in  the  two-dimensional  case  where  there  seem  to  be 
no  competing  algorithms. 

The  thinning  algorithm  can  also  be  implemented  more 
efficiently  at  the  cost  of  programming  complexity  and  by  using  a 
nonhomogeneous  bounding  process.  In  particular  the  method  can  be 
used  in  conjunction  with  the  special  algorithms  given  by  Lewis 
and  Shedler  (1976a,  1977). 

It  is  also  possible  to  extend  tha  method  of  thinning  to 
simulation  of  doubly  stochastic  or  conditioned  Poisson  processes. 

This  will  be  discussed  elsewhere. 
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